Wyzwanie konkursowe polega na predykcji wskaźnika bezrobocia na lata 2022-2023, na podstawie danych z 2000-2021. Analizę będziemy przeprowadzać przy użyciu języka R w programie RStudio, do obliczeń i przewidywań użyjemy m.in. z bibliotek takich jak stats, tseries, forecast, prophet, do wizualizacji - ggplot2, a do utworzenia ostatecznego raportu posłużymy się oprogramowaniem RMarkdown.
Wejściowe dane przyjmują formę szeregu czasowego - realizacji procesu stochastycznego, którego dziedziną jest czas - w tym przypadku po 12 miesięcy z 22 lat. Procesem stochastycznym \((X_t)_{t \in T}\) nazywamy rodzinę zmiennych losowych z pewnej przestrzeni probabilistycznej, przyjmującą wartości z przestrzeni mierzalnej.
Dane wskaźnika bezrobocia w latach 2000-2021 przedstawiają się w następujący sposób:
Z podziałem na lata:
Od razu zauważamy, że dane podlegają dostrzegalnemu, nieliniowemu trendowi - wartość wskaźnika spada na początku roku, w okolicy połowy roku wzrasta, by potem delikatnie spadać, lub utrzymywać się aż do października, a na koniec roku znowu wzrasta. Występuje oczywiście sezonowość o okresie 12 miesięcy
By sprawdzić podejrzenia wynikające z wizualnej obserwacji powyższych wykresów, przeprowadzamy dekompozycję STL (Seasonal Decomposition of Time Series by Loess). Ma ona wskazać składniki szeregu czasowego - jego trend, sezonowość oraz reszty.
library(forecast)
ts_stl <- ts(df$Wskaznik, frequency = 12, start = c(2000,12))
autoplot(mstl(ts_stl),col=ts_stl)
Powyższa metoda opiera się na przedstawieniu punktów szeregu czasowego (\(y_i, i \in T\)) jako suma komponentów sezonowości \(s_i\), trendu \(t_i\) i reszty \(r_i\): \[y_i = s_i + t_i + r_i\] oraz estymacja owych komponentów. [1]
Obserwując surowe dane, widzimy pewną anomalię w roku 2020 - wyraźny skok wartości wskaźnika bezrobocia w kwietniu w wyniku wybuchu pandemii COVID-19. To zaburzenie w danych w większości przypadków obniży jakość prognozy, ponieważ trend w tym okresie zostaje naruszony. Z tym problemem możemy poradzić sobie na kilka sposobów. Istnieje opcja podstawienia średniej wartości wskaźnika z każdego miesiąca do odpowiednich miesięcy z 2020. Można wziąć średnie globalne, lub jedynie z kilku ostatnich lat. Nie ma również większych przeszkód, by zupełnie pominąć ten rok w obliczeniach. Zbadamy także ideę, by użyć danych z lat 2000-2019, by “przewidzieć” wartości z 2020, a następnie dokonywać obliczeń przy użyciu nowych danych z 2020 do predykcji 2022-2023. Dokonamy analizy tych metod w rozdziale 3.
Jako pierwszą metodę predykcji wybraliśmy regresję liniową ze względu na miesiące. Wartości wskaźnika z n-tego miesiąca z lat 2000-2021 wyznaczają przewidywaną wartość wskaźnika z n-tego miesiąca na lata 2022 i 2023. Ze względu na podatność regresji liniowej na obserwacje odstające zdecydowaliśmy się usunąć rok 2020, ponieważ metoda ta nie wymaga ciągłości danych. [CITATION2]
library(stats)
library(tseries)
library(forecast)
predict_monthly <- function(data, month) {
monthly_data <- data %>% filter(M.c == month)
model <- lm(Wskaznik ~ Rok, data = monthly_data)
future_years <- data.frame(Rok = c(2022, 2023), M.c = month)
predictionsWith2020 <- predict(model, newdata = future_years)
return(data.frame(Rok = future_years$Rok, M.c = future_years$M.c, Wskaznik = predictionsWith2020))
}
predictionsWithout2020 <- lapply(1:12, function(m) predict_monthly(data, m))
predictionsWithout2020_df <- do.call(rbind, predictionsWithout2020)
combined_data <- bind_rows(data, predictionsWithout2020_df)
Ponieważ regresja liniowa jest ogólnym modelem, który nie jest dedykowany dla szeregów czasowych, zdecydowaliśmy się rozszerzyć te metodę o model ARIMA (Auto Regressive Integrated Moving Average). [CITATION3]
W modelu tym, część autoregresyjna (AR) jest zasadniczo formą regresji liniowej, gdzie bieżąca wartość szeregu czasowego jest modelowana jako liniowa kombinacja jego przeszłych wartości.
Ponieważ ARIMA również jest podatna na obserwacje odstające oraz wymaga ciągłości danych, zdecydowaliśmy się zastąpić rok 2020 predykcją modelu ARIMA na podstawie lat 2000-2019 i wykorzystać te dane aby przewidzieć wartości wskaźnika na podstawie lat 2000-2021 ze sztucznymi danymi z roku 2020. Mając swiadomość możliwości wystąpienia efektu kaskadowania błędów, porównamy to rozwiązanie z klasycznym zastąpieniem roku 2020 średnimi z pozostałych lat.
Tak wyglądają wartości z 2020, których użyjemy w predykcji lat 2022-2023:
fit_to_2019 <- auto.arima(ts_data)
forecast_2020 <- forecast(fit_to_2019, h=12)
Korzystając z funkcji adf.test sprawdzamy, czy spełniona jest stacjonarność modelowanych danych, która jest kluczowym założeniem modelu ARIMA.
adf.test(ts_data_doubleARIMA)
##
## Augmented Dickey-Fuller Test
##
## data: ts_data_doubleARIMA
## Dickey-Fuller = -9.5335, Lag order = 6, p-value = 0.01
## alternative hypothesis: stationary
adf.test(ts_data_mean)
##
## Augmented Dickey-Fuller Test
##
## data: ts_data_mean
## Dickey-Fuller = -9.4633, Lag order = 6, p-value = 0.01
## alternative hypothesis: stationary
KOMENTARZ JAKIE WYNIKI TESTOW
Zbadajmy następnie błędy ?????? co na skladowych wykresu
fit_doubleARIMA <- auto.arima(ts_data_doubleARIMA)
forecast_values_doubleARIMA <- forecast(fit_doubleARIMA, h=24)
fit_mean <- auto.arima(ts_data_mean)
forecast_values_mean <- forecast(fit_mean, h=24)
checkresiduals(fit_mean)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(1,0,1)(0,1,2)[12]
## Q* = 27.002, df = 20, p-value = 0.1352
##
## Model df: 4. Total lags used: 24
checkresiduals(fit_doubleARIMA)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(1,0,1)(0,1,2)[12]
## Q* = 27.762, df = 20, p-value = 0.1152
##
## Model df: 4. Total lags used: 24
Na powyższych wykresach błędów modelu widać, że zastosowanie podwójnej predykcji (drugi przypadek) jest teoretycznie lepsze niż model biorący rok 2020 jako średnią z pozostałych lat DLACZEGO. Jednakze, błędy w 2020 są praktycznie zerowe - wynika to z faktu, że rok ten jest dopasowany przez ten sam model. Mimo to zdecydowaliśmy się przyjąć wyniki z modelu korzystającego z podwójnej predykcji. LJUNG BOX TEST??
forecast_values_doubleARIMA <- data.frame(
Date = seq(as.Date("2022-01-01"), by = "month", length.out = 24),
Point_Forecast = forecast_values_doubleARIMA$mean[1:24]
)
forecast_values_mean <- data.frame(
Date = seq(as.Date("2022-01-01"), by = "month", length.out = 24),
Point_Forecast = forecast_values_mean$mean[1:24]
)
df_plot$Type <- "Actual"
forecast_values_doubleARIMA$Type <- "Double ARIMA"
names(forecast_values_doubleARIMA)[2] <- "Wskaznik"
forecast_values_mean$Type <- "Mean"
names(forecast_values_mean)[2] <- "Wskaznik"
combined_df <- rbind(df_plot[, c("Date", "Wskaznik", "Type")],
forecast_values_doubleARIMA,
forecast_values_mean)
Ostatecznie, dane przewidziane przy użyciu omówionych metod modeli regresyjnych przedstawiają się następująco:
ARIMA <- combined_df[combined_df$Type == 'Double ARIMA',]
Kolejną metodą przewidywania szeregu czasowego jest Prophet, przedstawiony przez Facebooka.[4] Na wyjściowym szeregu dokonujemy dekompozycji w następujący sposób: \[y(t) = g(t) + s(t) + h(t) + e_t\] W tym przypadku \(g(t)\) jest funkcją trendu reprezentującą nieokresowe zmiany wartości szeregu czasowego, \(s(t)\) jest funkcją zmian okresowych (np. miesięcznych), a \(h(t)\) to efekt świąt, występujących w nieregularnych odstępach czasu. Zakładamy, że błąd \(e_t\) ma rozkład normalny. funkcje \(g(t)\), \(s(t)\), \(h(t)\) są estymowane przy użyciu m. in. logistycznego modelu wzrostu oraz szeregów Fouriera, ściślej opisane w [4].
Podstawiając nasze dane uzyskujemy następujące wyniki predykcji:
library(prophet)
df_prophet <- df
df_prophet$date <- as.Date(paste(df$Rok, df$M.c, "01", sep = "-"))
df_prophet <- subset(df_prophet, select = c(date, Wskaznik))
colnames(df_prophet) <- c('ds', 'y')
model_prophet <- prophet(df_prophet)
future <- make_future_dataframe(model_prophet,
periods = 24,
freq = 'month')
forecast <- predict(model_prophet, future)
forecast[c('ds', 'yhat')] %>%
tail(n = 24)
## ds yhat
## 265 2022-01-01 105.59152
## 266 2022-02-01 100.91219
## 267 2022-03-01 98.81112
## 268 2022-04-01 97.07026
## 269 2022-05-01 97.14200
## 270 2022-06-01 97.91921
## 271 2022-07-01 99.66101
## 272 2022-08-01 99.95366
## 273 2022-09-01 99.73869
## 274 2022-10-01 99.50115
## 275 2022-11-01 101.79561
## 276 2022-12-01 103.00231
## 277 2023-01-01 105.79901
## 278 2023-02-01 101.15992
## 279 2023-03-01 99.35893
## 280 2023-04-01 97.09987
## 281 2023-05-01 97.20263
## 282 2023-06-01 98.16253
## 283 2023-07-01 100.04565
## 284 2023-08-01 100.28066
## 285 2023-09-01 99.94061
## 286 2023-10-01 99.66011
## 287 2023-11-01 101.96982
## 288 2023-12-01 103.25350
prophet_plot_components(model_prophet, forecast)
plot(model_prophet, forecast)
Na pierwszym wykresie widzimy krzywą trendu wyznaczoną przez model, łącznie z przewidzianymi ostatnimi latami wraz z przedziałem ufności (niebieskie pole na końcu krzywej), na drugim rysunku - uśrednione, ogólne roczne zmiany wskaźnika. Ostatni wykres przedstawia dodatkowo porównanie rzeczywistych danych (czarne punkty) z tymi estymowanymi przez Prophet (niebieska linia, wraz z przedziałem ufności).
Zbadajmy teraz jakość predykcji, kiedy pominiemy rok 2020 w rozważaniach:
forecast2[c('ds', 'yhat')] %>%
tail(n = 24)
## ds yhat
## 253 2022-01-01 104.48811
## 254 2022-02-01 99.89273
## 255 2022-03-01 97.55734
## 256 2022-04-01 95.63517
## 257 2022-05-01 95.70601
## 258 2022-06-01 96.53879
## 259 2022-07-01 98.34170
## 260 2022-08-01 98.68818
## 261 2022-09-01 98.51341
## 262 2022-10-01 98.31010
## 263 2022-11-01 100.63690
## 264 2022-12-01 101.82312
## 265 2023-01-01 104.36696
## 266 2023-02-01 99.81617
## 267 2023-03-01 97.72722
## 268 2023-04-01 95.69684
## 269 2023-05-01 95.67828
## 270 2023-06-01 96.49088
## 271 2023-07-01 98.35158
## 272 2023-08-01 98.70011
## 273 2023-09-01 98.48975
## 274 2023-10-01 98.33266
## 275 2023-11-01 100.65959
## 276 2023-12-01 101.92684
prophet_plot_components(model_prophet2, forecast2)
Obserwujemy, że po wyrzuceniu 2020, linia trendu zdecydowanie się wypłaszcza na koniec badanego okresu. Przez to, że model bierze pod uwagę współczynnik świąt, czyli nieregularnych skoków badanej zmiennej, w tym przypadku warto zostawić dane z 2020, zatem ostatecznie bierzemy pierwotną wersję prognozy.
PROPHET <- forecast['yhat'] %>%
tail(n = 24)
Metoda TBATS należy do popularnej grupy modeli statystycznych - modeli wygładzania wykładniczego. Do tej samej rodziny należy także STL, którą używaliścy do dekompozycji naszych danych.
Najczęściej używany model sezonowy przedstawia się w następujący sposób: [5]
\[y_t = l_{t-1} + b_{t-1} + s_t + d_t\] \[l_t = l_{t-1} + b_{t-1} + \alpha d_t\] \[b_t = b_{t-1} + \beta d_t\] \[s_t = s_{t-m} + \gamma d_t\] gdzie \(m\) to okres cykli sezonowych, \(d_t\) reprezentuje losowy szum w danych, \(l_t\), \(b_t\), \(s_t\) to komponenty odpowiednio: poziomu, trendu i sezonowości szeregu czasowego. Wartości \(\alpha\), \(\beta\) i \(\gamma\) są tak zwanymi parametrami wygładzającymi, a \(l_0\), \(b_0\), \(\{s_{1-m}, ..., s_0\}\) to zmienne wyjściowe. W tym modelu estymuje się właśnie zmienne wyjściowe oraz parametry. Warto dodać, że można estymować 2 składniki sezonowości \(s_t^{(1)}\), \(s_t^{(2)}\) wraz z parametrami \(\gamma_1\) i \(\gamma_2\), jednak w naszym przypadku bierzemy tylko jeden składnik - miesięczny.
Następnie dokonujemy modyfikacji tego modelu, przy uwzględnieniu transformacji Box-Cox, błędów modelu ARMA oraz T wzorców sezonowości. Na koniec, składniki sezonowości przedstawiamy jako ich trygonometryczną reprezentację opartą na szeregach Fouriera. Ze wszystkich składowych modelu powstała jego nazwa: T - trigonometrical, B - Box-Cox transformation, A - ARMA errors, TS - T seasonal patterns.
W celu ewaluacji jakości prognozy, będziemy porównywać ostatnie dane 2 lata z przewidzianymi wartościami na następne 2 lata, korzystając z dwóch metryk:
Po podstawieniu naszych danych otrzymujemy następujące wyniki:
ts_tbats <- msts(df$Wskaznik, seasonal.periods = 12)
model.ts_tbats <- tbats(ts_tbats)
model.ts_tbats.forecast <- forecast(model.ts_tbats, h = 24)
plot(model.ts_tbats.forecast, main = 'Prognoza TBATS',
ylab = 'Wskaznik')
model.ts_tbats.forecast$mean
## Jan Feb Mar Apr May Jun Jul
## 23 102.29247 98.84214 96.22711 95.65658 95.33452 96.74659 97.75088
## 24 104.10502 100.24764 97.32554 96.53332 96.03627 97.31860 98.21517
## Aug Sep Oct Nov Dec
## 23 98.91372 98.21906 98.66880 100.31739 102.23146
## 24 99.29118 98.52024 98.91193 100.51604 102.39415
W wykresie prognozy, niebieska linia oznacza przewidziane dane, wraz z przedziałami ufności na poziomie istotności \(0.2\) oraz \(0.05\). Przewidziane wartości zostały wypisane w tabeli pod nim.
df.test <- tail(df$Wskaznik, n = 12 * 2)
ts_tbats.predict <- predict(model.ts_tbats.forecast, df.test)
ts_tbats.test_vs_predicted <- data.frame(df.test, model.ts_tbats.forecast$mean)
matplot(ts_tbats.test_vs_predicted, type = 'l', lty = 1:2, col = 1:2, ylab = 'Przypadek nr 1')
Na tym rysunku widzimy porównanie wartości z ostatnich dwóch lat (czarna linia) z przewidzianymi wartościami na następne 2 lata (czerwona przerywana linia). Obserwujemy tutaj zauważalne niedopasowanie krzywych na początku okresu - jest to następstwo wzięcia pod uwagę odstających wartości z roku 2020.
MAPE.error(ts_tbats.test_vs_predicted$df.test,
ts_tbats.test_vs_predicted$model.ts_tbats.forecast.mean)
## [1] 2.085391
RMSE.error(ts_tbats.test_vs_predicted$df.test,
ts_tbats.test_vs_predicted$model.ts_tbats.forecast.mean)
## [1] 6.855342
Ostatecznie, liczymy wartości błędów MAPE i RMSE - są one zawyżone z wyżej wymienionego powodu.
Aby zmiejszyć wartości błędów, dokonujemy analizy bez 2020:
## MAPE: 1.168697
## RMSE: 3.977294
W przypadku wyrzucenia danych z 2020, zgodnie z oczekiwaniami, widzimy znaczny spadek błędów MAPE i RMSE. Zwróćmy uwagę, że do ewaluacji prognozy wzięliśmy lata 2019 i 2021 zamiast 2020-2021. Sprawdzimy jeszcze predykcję dla przypadków, gdy zamiast wartości z tego roku podstawimy średnią globalną wartość wskaźnika:
## MAPE: 1.009071
## RMSE: 2.643147
Po raz kolejny, wartości błędu zmiejszyły się. Na koniec zbadajmy zachowanie błędów, gdy zamiast 2020 weźmiemy średni wskaźnik z ostatnich 5 lat:
## MAPE: 1.06453
## RMSE: 3.255959
Najmniejsze wartości błędów otrzymujemy w przypadku podstawienia globalnej średniej za wartości z odstającego roku. Decydujemy jednak uwzględnić ostatni przypadek w ostatecznej predykcji, gdyż nie chcemy brać takich wyników, które będą zbyt zbliżone do wartości z ostatnich dwóch lat, ponieważ po raz kolejny, trend byłby zbyt spłaszczony na końcu okresu. Ostatecznie:
TBATS <- model.ts_tbats4.forecast$mean %>% as.data.frame()
Wszystkie wyniki wygenerowane przez omówione modele predykcji przedstawiają się w następujący sposób:
## Date ARIMA Prophet TBATS
## 1 2022-01-01 102.82608 105.59152 102.15098
## 2 2022-02-01 98.62021 100.91219 98.84227
## 3 2022-03-01 96.13074 98.81112 96.15523
## 4 2022-04-01 95.42150 97.07026 95.20574
## 5 2022-05-01 95.86357 97.14200 94.80567
## 6 2022-06-01 95.77652 97.91921 96.32381
## 7 2022-07-01 97.74766 99.66101 97.41151
## 8 2022-08-01 98.41741 99.95366 98.46648
## 9 2022-09-01 97.32115 99.73869 97.71175
## 10 2022-10-01 97.46430 99.50115 98.19079
## 11 2022-11-01 99.36107 101.79561 99.75118
## 12 2022-12-01 100.32231 103.00231 101.52969
## 13 2023-01-01 103.52965 105.79901 103.32886
## 14 2023-02-01 99.01594 101.15992 99.75306
## 15 2023-03-01 96.44852 99.35893 96.86344
## 16 2023-04-01 95.57216 97.09987 95.76631
## 17 2023-05-01 96.07756 97.20263 95.25199
## 18 2023-06-01 96.01987 98.16253 96.68639
## 19 2023-07-01 98.06441 100.04565 97.70474
## 20 2023-08-01 98.71224 100.28066 98.70353
## 21 2023-09-01 97.60887 99.94061 97.89990
## 22 2023-10-01 97.69950 99.66011 98.34201
## 23 2023-11-01 99.70443 101.96982 99.87406
## 24 2023-12-01 100.64744 103.25350 101.62973
Ostateczne wyniki otrzymujemy biorąc średnią dla każdego miesiąca ze wszystkich metod:
Wyznaczona przez nas predykcja, łącząca w sobie modele ARIMA, Prophet oraz TBATS jako średnia arytmetyczna wszystkich prognoz prezentuje się następująco:
## Date Wskaznik
## 1 2022-01 103.52286
## 2 2022-02 99.45822
## 3 2022-03 97.03237
## 4 2022-04 95.89917
## 5 2022-05 95.93708
## 6 2022-06 96.67318
## 7 2022-07 98.27339
## 8 2022-08 98.94585
## 9 2022-09 98.25720
## 10 2022-10 98.38541
## 11 2022-11 100.30262
## 12 2022-12 101.61810
## 13 2023-01 104.21917
## 14 2023-02 99.97631
## 15 2023-03 97.55696
## 16 2023-04 96.14611
## 17 2023-05 96.17739
## 18 2023-06 96.95626
## 19 2023-07 98.60494
## 20 2023-08 99.23214
## 21 2023-09 98.48312
## 22 2023-10 98.56721
## 23 2023-11 100.51611
## 24 2023-12 101.84355
Cleveland, R. B., Cleveland, W. S., McRae, J. E., & Terpenning, I. J. (1990). STL: A seasonal-trend decomposition procedure based on loess. Journal of Official Statistics, 6(1), 3–33.
regresja
arima
Taylor SJ, Letham B. Forecasting at Scale. The American Statistician. 2018
De Livera A.M., Hyndman R.J., Snyder R.D., Forecasting time series with complex seasonal patterns using exponential smoothing. Journal of the American Statistical Association 2011, 106, 1513–1527.